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The reduced density matrix of an interacting system can be used as the basis for a truncation scheme, or in an 
unbiased method to discover the strongest kind of correlation in the ground state. In this paper, we investigate 
the structure of the many-body fermion density matrix of a small cluster in a square lattice. The cluster density 
matrix is evaluated numerically over a set of finite systems, subject to non-square periodic boundary conditions 
given by the lattice vectors R| = (Ri x ,Ri y ) and R2 = (R2x,R2y>- We then approximate the infinite-system 
cluster density-matrix spectrum, by averaging the finite-system cluster density matrix (i) over degeneracies in 
the ground state, and orientations of the system relative to the cluster, to ensure it has the proper point-group 
symmetry; and (ii) over various twist boundary conditions to reduce finite size effects. We then compare the 
eigenvalue structure of the averaged cluster density matrix for noninteracting and strongly-interacting spinless 
fermions, as a function of the filling fraction h, and discuss whether it can be approximated as being built up 
from a truncated set of single-particle operators. 



I. INTRODUCTION 



The density matrix is a very useful tool in the numer- 
ical study of interacting systems. Besides being used in 
the Density-Matrix Renormalization Group (DMRG) 1 and 
its higher-dimensional generalizations, 2 the density matrix is 
also used as a diagnostic tool in the Contractor Renormaliza- 
tion (CORE) method for numerical renormalization group in 
two dimensions, 3 and forms the basis of a method to identify 
the order parameter related to a quasi-degeneracy of ground 
states. 4 

In previous work, 5 we extended the results of Chung and 
Peschel 6 to write the density matrix (DM) of a cluster of Nc 
sites cut out from a system of noninteracting spinless fermions 
in d dimensions as the exponential of a quadratic operator, 
called the pseudo-Hamiltonian, as it resembles the Hamilto- 
nian of a noninteracting system. That result was then applied 
in numerical studies of noninteracting spinless fermions in 
one dimension, to better understand how the distribution of 
cluster DM eigenvalues scale with Nc, and to explore the pos- 
sibility of designing truncation schemes based on the pseudo- 
Hamiltonian. 7 We believe truncation schemes such as that de- 
scribed in Ref.7 will be helpful to the choice of basis states in 
renormalization groups such as CORE. 

Thus, some questions motivating the present paper were: 
(i) does the density matrix of an interacting Fermi-liquid sys- 
tem resembles that of a noninteracting one? (ii) can we apply 
our exact result in Ref.5 to two dimensions as well as for one 
dimension? (iii) is it numerically practical to compute this 
sort of density matrix in a fermion system. To answer these 
questions, we investigated a spinless analog of the extended 



Hubbard model, given by the Hamiltonian 

//= - f Z c / c v + y Z w ' n > (L1) 

(i.j) (i,f> 

in the limit of V — > 00, so that fermions are not allowed to 
be nearest neighbors of each other. This model is chosen for 
two reasons: (i) for a given number of particles, the V — > 00 
Hilbert space is significantly smaller than the V < 00 Hilbert 
space, and we can work numerically with larger system sizes; 
and (ii) the model, in spite of its simplicity, has a rich zero- 
temperature phase diagram, 8-10 where we find practically free 
fermions in the limit n <k 1, and an inert solid at half-filling 
h — i. As the filling fraction approaches quarter- filling from 
below, h — > j , the system becomes congested, highly cor- 
related, but is nonetheless a Fermi liquid, perhaps with ad- 
ditional orders that are not clear in small systems. Slightly 
above quarter-filling, the dense fluid and inert solid coexists, 
while slightly below half-filling, the system is expected to sup- 
port stable arrays of stripes. 

To probe this rich variety of structures in the ground state at 
different filling fraction «, we describe in Section II how the 
reduced DM of a small cluster, with the appropriate symme- 
try properties, can be calculated from a finite non-square sys- 
tem subject to twist boundary conditions. Then in Section III, 
we investigate in great details the cluster DM spectra of the 
noninteracting system, particularly on how to handle finite 
size effects in the numerics, for comparison with the cluster 
DM spectra of a strongly-interacting system, presented in Sec- 
tion IV. Finally, in Section V, we summarized our findings, 
and discuss the prospects of designing an Operator-Based DM 
Truncation Scheme for interacting systems, at some, if not at 
all, filling fractions. 
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II. FORMULATION 

In this section, we give the theoretical formulations and 
describe the numerical tools needed to investigate the clus- 
ter DM spectra of noninteracting and strongly-interacting sys- 
tems of spinless fermions in two dimensions. In Section II A, 
we give the matrix elements of the DM of a small cluster em- 
bedded within a larger, but still finite, system. These matrix 
elements are obtained by tracing out degrees of freedom ex- 
ternal to the cluster, starting from the ground-state wave func- 
tion of the system, obtained through exact diagonalization. In 
Section II B, we describe how our finite systems can be de- 
fined with nonsquare periodic boundary conditions, and how 
we make use of the translational invariance of both noninter- 
acting and strongly-interacting models to reduce the computa- 
tional efforts in exact diagonalization. In Section II C, we de- 
scribe several averaging apparatus required to obtain a handle 
on the infinite-system spectra of the cluster DM, and then in 
Section II D, we describe a classification scheme for the one- 
particle and multi-particle eigenstates of the cluster DM that 
makes the symmetry of the underlying square lattice explicit. 

A. Cluster Density Matrix 

The DM pc of a cluster cut out from a larger system is a 
density operator which gives the expectation 

CflApP) = (A) = Tr c p c ^ (2.1) 

for any observable A local to the cluster, when the larger sys- 
tem is in its ground state The cluster DM pc can be cal- 
culated from the ground-state DM 

p = mm (2.2) 

of the system, by tracing out degrees of freedom outside of the 
cluster. We write this as 

p c =Tr E p, (2.3) 

where the subscript E denotes a trace over environmental de- 
grees of freedom. 

Since a cluster is a collection of sites identified in real 
space, it is natural to choose as a many-body basis the real- 
space configurations. For a finite two-dimensional system 
with N sites, we label the sites 7 = 1 through j - N, so that 
for any pair of sites (xj^y^) and (xj 2 ,yj 2 ), we have xj t < Xj 2 
and y^ < yj 2 if ji < ]2- We then distinguish between sites 
within the cluster, of which there are Nr of them, (x*:,y,c), 
(x.c, v -c), . . . , (xf ,y f ), and sites outside of the cluster, of 

which there are Ne = N - Nc of them, (.*,■£, y^), (x^y^), 
. . . , (x :e , y :e ). We think of the Nf sites outside the cluster as 

Jn e j Jn e 

constituting the environment to the cluster. 

We work with the configuration basis states |j) = 
\jih ' ' ' jp)> where j\ < ■ ■■ < jp are the P occupied sites 
in the system. These can be thought of as a direct product of 
the configuration basis states of the cluster |1) = ■ ■ ■ lp c ), 



where l\ < ■ ■ ■ < lp c are the Pc occupied sites within the clus- 
ter, and the configuration basis states of the environment |m) = 
\m\ni2 ■ ■ ■ nip E ), where m\ < ■■■ < nip E are the Pe — P - Pc 
occupied sites in the environment. Here, we have the occupied 
sites of the system ., jp} = . . . ,lp c ] U {mi, . . ., mp E ) 
being the union of the occupied sites in the cluster and in the 
environment, with the site indices I and m resorted in ascend- 
ing order to give the site indices j. 

In terms of the configuration basis of the system, the 
ground-state wave function of the system can be written as 

\^) = ^^) = Y J ^c) i ---c) p \0), (2.4) 
j j 

where is the amplitude associated with configuration |j), 
and Cj, cj are fermion annihilation and creation operators act- 
ing on the site {Xj,yj). We can also write the ground-state 
wave function as 

^> = XZ ( - 1)/<j;l ' m,>p, ' m|1>|in) 

1 m 

= ZZ ( - 1)/ti;lm)>I,| ' mX (2 - 5) 

1 m 

c mi ' ' ' c Lp E c i, " ' c J Pc 10) ' 

in terms of the direct product of configuration bases of the 
cluster and the environment, where c t and are fermion an- 
nihilation and creation operators acting on site (xi,yi) within 
the cluster, and c m and c] n are fermion annihilation and cre- 
ation operators acting on site (x m ,y m ) within the environment. 
In (2.5), the amplitude m = is taken directly from the ex- 
pansion in (2.4), while the factor (— 1)/0A»0 accounts for the 
fermion sign incurred when we reorder the operator product 
4, ' ' ' c m PE c { ■ ■ ■ c l c to § et the operator product ■ ■ ■ c] ? . 
Similarly, the ground-state DM in (2.2) can be written as 

/'-XXw; <>><«><■,...<> (2.6) 
j j' 

using the system-wide configuration basis, or as 

l.m I',m' 

c m"- c 'LA''' C U mX 

<0|c„ ■■■c„c , - c (2.7) 

r c r E 

using the direct-product basis between cluster configurations 
and environment configurations. 

Performing the trace over the environment as prescribed in 
(2.3), we find the fermion cluster DM pc to be 

p c = 2Z ( " 1)/ti;1 ' m)+/(i ' ;1 ' m)x 

1,1' m,m' 

<•; •••<•; ioxoi c, ■■■c r (2.8) 

>1 'P C lp,^ l! 
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Its matrix elements are 

<\\Pc\V) = ^Y J (-l) fQXm)+fQ '*' m ' ) x 

m m 

Tlm^W- (2-9) 

These matrix elements can be computed naively by run- 
ning over all possible pairs of cluster states |1) and |1'), and 
performing the sums over m and m' as they appear in (2.9), 
looking up the fermion-sign factors (-l)/<J;l,ni) anc j amplitudes 

m as and when they are needed. We call this the naive al- 
gorithm. Alternatively, we can also reorganize the fermion- 
sign-factor-adjusted amplitudes (— ly^'^'Wy^ into a matrix 
*P, whose rows are associated with the cluster configurations 
[1), and whose columns are associated with the environment 
configurations |m), the matrix pc can be computed directly by 
matrix multiplication as 

p c = ** f , (2.10) 

representing a collection of inner products. We call this the 
presorted inner-product algorithm. We compare and analyze 
the computational complexity of these two algorithms in Ap- 
pendix A. In the numerical studies presented in Section IV, 
we use the pre-sorted inner-product algorithm exclusively. 

From (2.10), we see that the cluster DM pc is manifestly 
hermitian, and thus all its eigenvalues {w} are real (and in 
fact nonnegative). When obtained from the wave function of 
a state with definite particle number, pc has no matrix ele- 
ments between cluster states containing different number of 
particles, and thus the eigenstates \w) of pc can be organized 
into sectors, corresponding to Pc - 0, 1, ... , Pc.max particles 
within the cluster. For the rest of this paper, we would re- 
fer to the eigenvalues of pc generically as its weights, since 
these have a natural probabilistic interpretation. A Pc-particle 
weight of pc is therefore an eigenvalue corresponding to an 
eigenstate containing Pc particles within the cluster. 

B. System Definition and Translational Invariance 

For noninteracting spinless fermions on an infinite square 
lattice, it is possible to compute the cluster DM pc starting 
from the Fermi sea ground state, through the evaluation and 
diagonalization of Gc- For an interacting system, we need 
to compute pc starting from p in (2.2), the latter we obtain 
through exact diagonalization on a finite system. We define 
the a finite system relative to the infinite square lattice in terms 
of the lattice vectors Ri and R2, as shown in Fig. 1, such that 
N - z ■ (Ri x R 2 ) = RixRiv - RixRiv > is the number of 
lattice sites within the system. 

If we impose periodic boundary condition such that R + 
mRi + 11R2 = R, then in the exact diagonalization to obtain 
I*!*) we can take advantage of translational invariance through 
the use of the Bloch states 

Lj ,k>= -^yy ik - R r R |j () > (2.H) 




FIG. 1: Definition of system to be exactly diagonalized in terms of 
the lattice vectors R[ and R 2 . We shall denote such a system as 
R[ xRj. In this example shown, the system is (5, 1) x (1,4). 

as our computational basis. 9 In this Bloch state, the config- 
urations {Tr jo)} are all related to the generating configura- 
tion jo) by the lattice translations Tr associated with dis- 
placement R, while k are wave vectors allowed by the bound- 
ary conditions. Any configuration within the collection of 
translationally-related configurations {Tr |j }} can serve as the 
generating configuration, but we pick the one with the least 
sum of indices of occupied sites. 

Working with finite non-square systems introduces several 
complications. First of all, we sometimes end up with degen- 
erate ground states which suffer from symmetry-breaking not 
found in the true infinite-system ground state. However, be- 
cause the point symmetry group of our non-square finite sys- 
tem is only a subgroup of the square lattice point symmetry 
group, the finite-system ground-state manifold is not invariant 
under all square lattice symmetry operations. Thirdly, when 
working with finite systems, we introduce systematic devia- 
tions which are collectively known as finite size effects. We 
identify the three primary sources of finite size effects as (i) 
finite domain effect, which has to do with the fact that the 
small set of discrete wave vectors allowed are not adequately 
representative of the continuous set of wave vectors on the in- 
finite square lattice; (ii) shell effect, which has to do with the 
fact that the set of discrete wave vectors allowed are organized 
by symmetry into shells in reciprocal space, each of which can 
be partially or fully filled in the many-body ground state; and 
(iii) shape effect, which has to do with the detailed shape of 
the non-square system we introduced. 



C. Averaging 

1. Degeneracy Averaging 

To eliminate these numerical artefacts, we adopted three av- 
eraging devices. First, we average over the Do-fold degenerate 
ground-state manifold. Our first motivation for doing so is as 
follows: if ^ is the point symmetry group of the square lattice, 
and its subgroup G is the point symmetry group of the Ri x R2 
system, then we will find that the cluster density matrices 

Pa = Tfg |\Pj) (T*,! , (2.12) 
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one for each wave function ^F, within the Do-fold degenerate 
ground-state manifold, are not invariant under G, much less 
We remove this artificial symmetry breaking by calculating 
the degeneracy-averaged cluster DM 



j d 



(2.13) 



over the cluster density matrices p c , within the ground-state 
manifold. This degeneracy-averaged cluster DM is invariant 
under G. 

A second motivation for such a mode of averaging over the 
degenerate ground-state manifold of the finite system is that 
thermodynamically, given the pure state density matrices pcj, 
with energy eigenvalue we typically construct the canoni- 
cal ensemble DM as 



(2.14) 



where Z(p) - Yji e ^ E ' i s me canonical partition function. 
States within a degenerate manifold have the same energy, 
and therefore contribute equally to the thermodynamic DM 
PciP)- In the limit of /3 — > oo, the usual thermodynamic ar- 
gument is that pure states decouple from one another, and we 
treat their respective density matrices independently, except 
for those states which are degenerate. Because they appear 
with the same Boltzmann weight whatever the inverse temper- 
ature is, we should still treat the uniform combination instead 
of the individual density matrices in the limit of /3 — > co, 



2. Orientation Averaging 

The second averaging device involves an average over the 
orientation of the finite non-square system relative to the 
underlying square lattice. This averaging restores the ( S- 
symmetry to the averaged ground state. In principle, this 
requires us to compute pc for a group of four systems: 
(R u ,Ri y ) x {R 2x ,R 2y \ (R2v,Ri x ) x (Ri y ,Ri x ), (Ru,-Ru) x 
(-Ri x ,R2y) and (R 2y , -R 2x ) x (-Ri y , R u ). 

However, if the cluster whose DM we are calculating is in- 
variant under the action of ^ , this averaging can be achieved 
by computing 



1 



Pc 



(2.15) 



where g e G is a point group transformation of the square 
lattice, U(, is the unitary transformation of the cluster Hilbert 
space associated with g, and D(^) is the order of . 



3. Twist Boundary Conditions Averaging 

After these two averagings, the cluster DM has the full sym- 
metry (including translations) of the underlying square lattice, 



but finite size effects remain. We eliminate these as much as 
we can 11 with the third averaging device, twist boundary con- 
ditions averaging. 12 The usual way to implement twist bound- 
ary conditions is to work in the boundary gauge, keeping the 
Hamiltonian unchanged, and demanding that 



R+R, 



-i0-Ri 



R+R, 



i^R, 



(2.16) 



where R[ and R2 are the lattice vectors defining our finite sys- 
tem, R = (R x , R y ) is a site within the system, and <f> = (<p x , 4>y) 
is the twist vector parametrizing the twist boundary condi- 
tions. 

In choice of gauge (2.16), the Hamiltonian (1.1) is not man- 
ifestly invariant under translations. However, we can continue 
to block-diagonalize it using the Bloch basis states defined in 
(2.11), provided the set of allowed wave vectors k are shifted 
relative to ko for the usual periodic boundary conditions by 
the twist vector <f>, i.e. 



k = kn + - 



(2.17) 



The other natural way to implement twist boundary condi- 
tions is in the bond gauge, where we make the substitution 



in the Hamiltonian, but demand that 



R+R, 



R- 



R+R, 



(2.18) 



(2.19) 



where R[ and R2 are the lattice vectors defining our finite sys- 
tem. Now the Hamiltonian (1.1) is manifestly invariant under 
translations, and we can bloch-diagonalize it using the Bloch 
basis states defined in (2.11), with the same set of allowed 
wave vectors k = ko as for the usual periodic boundary con- 
ditions. 

Exact diagonalization can be performed in any gauge, but 
we chose to do it in the bond gauge, because in this gauge, the 
Bloch basis states [jo, k) defined in (2. 1 1) can be used as is to 
block diagonalize the Hamiltonian at all twist vectors <f>. This 
gives us the ground-state wave function ^(bond)) in the bond 
gauge. In the boundary gauge, or any other gauges, appropri- 
ate gauge transformations must be applied to |jo, k) before we 
can use this Bloch basis to block diagonalize the Hamiltonian. 
Because of this, the computational cost for exact diagonaliza- 
tion incurred in the bond gauge is fractionally lower than in 
other gauges. 

We can also calculate the correlation functions 
(0(R)0'(R + AR)> (of which the cluster DM is a func- 
tion of) in any gauge, with appropriately-defined covariant 
observables O = UOU', where O are the 'physical' observ- 
ables we would use when there is no twist in the boundary 
conditions. In the boundary gauge, these covariant observ- 
ables 0-0 and O' - O' are particularly simple, except 
when the displacement vector AR crosses the boundaries of 
our system. For our purpose of calculating the cluster DM, 
this situation occurs only when the cluster itself straddle 
the system boundaries. Therefore, with the cluster properly 
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nested within the system, we chose to perform twist boundary 
conditions averaging in the boundary gauge. 

In the boundary gauge, the cluster DM is obtained by trac- 
ing down the ground-state wave function ^(boundary)). We 
can get this wave function from |*F(bond)) by applying the 
gauge transformation 



f 



|n) -> e -''2R»R R ^ |n), 



(2.20) 



where |n) is an occupation number basis state, with occupation 
«r on site R. 

Now, averaging over twist vectors is the same as integrat- 
ing over the Brillouin Zone, so we perform twist boundary 
conditions averaging over a uniform grid of Monkhorst-Pack 
points with order q. 13 For rectangular systems (L x , 0) x (0, L y ), 
the First Brillouin Zone is sampled by varying the two in- 
dependent twist angles between -n/L x < <p x < +7t/L x and 
-n/Ly < <p y < +n/L y . For non-square systems, the two inde- 
pendent twist angles <p\ and <p2 are defined by 



J<t>*2 



(2.21) 



The twist vector <p is then related to the independent twist an- 
gles —Tt < (p\ < +n and —n<<p2< +n by 



^ n 4- ^ n 

Lti In 



(2.22) 



where 



2tt 



Ql = —(Rir.-Rixl 



2tt 

Qi = — (-RiyM 



are the primitive reciprocal lattice vectors of our non-square 
system. For such systems, the First Brillouin Zone is a 
parallelogram on the (p x -<f> y plane, so the uniform grid of 
Monkhorst-Pack points are imposed on the square domain 
(— 7T, +n) x (—7t, +7t) on the <f>i-<j>2 plane instead. 



D. Classifying States of the Cluster 

The point symmetry group of the square lattice is the di- 
hedral group D4, which has eight elements. 14 The five irre- 
ducible representations of this group are A\, A 2 , B\, B2 and E. 
For the cross-shaped cluster shown in Fig. 2, there is a one-to- 
one correspondence between these five irreducible representa- 
tions and the one-particle states, but instead of labeling these 
one-particle states as |Aj), A2), |Z?i), I-S2) an d \E), we adopt 
an angular momentum-like notation, 



|A,) = ^=(1,1, 1,1,1), 



| 5 _) = |A 2 )= ^=(1,1,-1,1,1), 
\Px) = \B 1 ) = ^(1,0,0,0,-1), 
|p,> = l«2> = ^(0,1,0,-1,0), 
\d) = \E} = 1(1,-1,0,-1,1). 



that would make clear the structure of these one-particle 
states. We find it more convenient to work with the one- 
particle states 



l-v> 



1 



ff + > + 4=|j_> = i(i,i,o,i,i), 



V2 |J + / ' V2 

s) = (0,0,1,0,0) 
1 



■1), 



\ P+ ) = -fe\P X )+^\Py)=±(l,i,0,-i, 

\P-)=^\p X )-^\ P y) = \(l,-i,0,i,-V,. 



(2.25) 




FIG. 2: The five-site, cross-shaped cluster whose cluster DM we are 
interested in calculating for both a system of noninteracting, as well 
as strongly-interacting spinless fermions. 

For a cluster DM possessing the full point group symmetry 
of the square lattice, the one-particle state \d) is constrained 
by symmetry to always be an eigenstate of p c . We call the 
associated eigenvalue w t i the weight of \d). Furthermore, the 
one-particle states \p x ) and \p y ) are also equivalent under the 
square lattice symmetry, and hence their weights w Px and w p 
are equal. We call this doubly-degenerate one-particle weight 
(2.23) w p . On the other hand, the s one-particle eigenstates of pc are 
in general not \s+), |s_) or \s), \s), but some admixture of the 
form 



\si) — cos6*|i) + sin 6\s) , 
|i2) = - sin 6 \s) + cos 8\s) . 



(2.26) 



We call their corresponding weights w s , and w Sl respectively. 

We can then extend this angular-momentum-like notation to 
multi-particle states of the cluster. Though the quantum num- 
bers used to label the one-particle states are, strictly speaking, 
not angular momentum quantum numbers, we apply the rules 
of angular momenta addition as if they were to write down 
the angular-momentum-like quantum numbers for the multi- 
particle states. For example, for the two-particle states of the 
cluster, we have 



\S) = \p + p-), \S') = \ss), 
\P + ) = \sp + ), \P' + ) = \sp + ), \P'l) = \p_d), 
\P-) = \spJ), \P'_) = \sp-), \P'l) = \p + d), 
\D) = \sd), \D') = \sd). 



(2.27) 



(2.24) HI. NONINTERACTING SYSTEM 

In preparation for our main calculations on the strongly- 
interacting system, we investigated in great details the cluster 
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DM spectra of a system of noninteracting spinless fermions 
described by the Hamiltonian 

H = -tY u c]c j . (3.1) 

<u> 

An analytical formula for the cluster DM pc is known for this 
system, 5 using which we can obtain the spectrum of p c for any 
system size. We take advantage of this analytical formula to 
calculate pc for an infinite system of noninteracting spinless 
fermions. 

However, our goal in calculating the cluster DM spectrum 
of noninteracting spinless fermions is to compare it against the 
cluster DM spectrum of interacting spinless fermions. For the 
latter system, we can only calculate — sans approximations 
— the cluster DM from the exactly-diagonalized ground-state 
wave function of finite systems. To make this comparison be- 
tween noninteracting and interacting spinless fermions more 
meaningful, we compute their cluster DMs for the same se- 
ries of finite systems. In Section III A, we describe how the 
infinite-system and finite-system cluster DMs for noninteract- 
ing spinless fermions are calculated. 

In Section IIIB, we observe that the noninteracting finite- 
system cluster DM spectra are contaminated by various finite- 
size effects. These finite-size effects also arise when we com- 
pute the interacting finite-system cluster DM spectra, so we 
want to learn how to deal with them. Clearly, the effectiveness 
of various techniques in reducing finite-size effects can be best 
gauged by applying them to finite systems of noninteracting 
spinless fermions, since we will then be able to compare the 
results from the various techniques against the infinite-system 
limit. The simplest antidote to the various finite-size effects 
is to use a larger finite system. As expected, finite-size ef- 
fects do become less and less important as the size of the 
system is increased. Unfortunately, based on comparisons of 
the finite-system cluster DM spectra with the infinite-system 
cluster DM spectrum, we realized that we would need to go 
to system sizes of a few hundred sites in order for the finite- 
system cluster DM spectra to be decent approximations of the 
infinite-system cluster DM spectrum. 

Since such system sizes are not practical for exactly di- 
agonalizing the strongly-interacting system given by (1.1), 
we look into the method of twist boundary conditions aver- 
aging in Section IIIC. This method involves averaging the 
cluster DM spectra over various phase twists introduced into 
the periodic boundary conditions imposed on a given finite 
system. For noninteracting spinless fermions, we find that 
twist boundary conditions averaging reduces finite domain 
and shell effects in the cluster DM spectra, which then approx- 
imate the infinite-system cluster DM very well. As a matter of 
standardization, we apply the method of twist boundary con- 
ditions averaging to both noninteracting and interacting spin- 
less fermions, and compare their twist boundary conditions- 
averaged cluster DM spectra. 



A. Calculating the Noninteracting Cluster DM 

Instead of the general formalism presented in Section II A, 
for the system of noninteracting spinless fermions we calcu- 
late the cluster DM weights using the exact formula 

p c = det(l - G c )exp{^ log e [G c (l - G c )] y cjc.} (3.2) 

ij 

obtained in Ref. 5, which relates the DM pc of a cluster of 
sites and the cluster Green-function matrix Gq- The matrix 
elements of Gq are given by 

G C (R,R') = 0F|4c R ,|¥> = i Yj e lk ' (R - R,) , (3.3) 

k filled 

where \W) is the ground state of the system, and R, R' are 
sites within the cluster. The corollary of (3.2) is that, if Ai 
is an eigenvalue of the cluster Green-function matrix Gc, the 
corresponding one-particle weight of pc is 

Wl = Ai Y\(l ~ Ar). (3.4) 

To calculate the infinite-system spectra of pc, we convert 
the sum over k„ in (3.3) into an integral 

G C (R-R')= f ^ lk ' (R - R ' , (3.5) 

over those wave vectors k bounded by the Fermi surface 
e(k) = €f, where €f is the Fermi energy. On an infinite 
square lattice with unit lattice constant, the dispersion relation 
is given by 

e(k) = -2 (cos k x + cos k y ) . (3.6) 

We then obtain the infinite-system cluster Green-function ma- 
trix eigenvalues A Si (cf), A p (€f), A^iep) and A s ,(€f) as func- 
tions of the Fermi energy ep by numerically integrating (3.5), 
and diagonalizing the resulting infinite-system cluster Green- 
function matrix Gc(ff). For the same set of Fermi energies, 
we also integrate 

over the Fermi surfaces to find the corresponding filling frac- 
tions. The one-particle infinite-system cluster DM weights 
w Si (€f), w p (6f), w^(ep), and w J2 (ep) are then calculated using 
(3.4). 

To calculate the cluster DM spectra for a finite system of N 
sites with P noninteracting particles, we determine the set of 
wave vectors {k„}^ =1 with the lowest single-particle energies 

£k„ = -2(cos£„ iA . + cos (3.8) 

We then evaluate the finite-system cluster Green-function ma- 
trix elements in (3.3) by summing over these occupied wave 
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vectors {k„}^ =1 . Following this, we diagonalize the finite- 
system cluster Green-function matrix Gc(P) to obtain the 
eigenvalues A Sl (P), A p (P), Ad(P), and A Sl (P), and therefrom 
the one-particle finite-system cluster DM weights w Sl (P), 
w p (P), Wd(P), and w s ^(P) using (3.4). By varying the number 
P of noninteracting particles in the finite A^-site system, we 
determined the finite-system cluster DM one-particle spectra 
for the filling fractions n = P/N accessible to the finite system. 



B. Finite Size Effects and the Infinite-System Limit 

Imposing the usual periodic boundary conditions, we calcu- 
lated the finite-system spectra of pc for several small systems 
ranging from N = 13 to N = 20 sites. For these small systems 
sizes, we find that it is impossible to say anything meaningful 
about the dependence on the filling fraction h = P/N for the 
cluster DM weights because of the finite size effects. Using 
the relation (3.2) for noninteracting systems, we investigated 
the effect of system size on the spectrum of pc for a series of 
system (4p, —p) x (p, 4p), 1 < p < 8, with the same shape. As 
the system size is increased, we find that the cluster DM spec- 
trum approaches the infinite-system limit, as shown in Figure 
3. 

For this series of systems, the infinite-system limit is more 
or less reached at around p = 4 (272 sites), based on compar- 
ison with the infinite-system limit itself. We can also arrive at 
this estimate by looking at the convergence of the one-particle 
weights alone. More importantly, we find that the shell effect 
affects weights of different symmetry differently: w Sl is al- 
most unaffected, while wj is the most severely affected. Shell 
effect persists in even up to a system size of 1088 sites (for 
P = 8). 

We also looked at w Sl (h), which is almost unaffected by 
shell effect, for several systems with between 200 to 300 
sites of different shapes. For systems of these sizes, the fi- 
nite domain effect is negligible, but we find w Sl (h) from fi- 
nite systems of different shapes differing very slightly from 
the infinite-system limit, and also from each other. Since we 
expect systems of different shapes to all approach the same 
infinite-system limit, we attribute these very small deviations 
to the shape effect. Based on more extensive numerical studies 
(see Chapter 4 of Ref. 1 1) not reported in this paper, we know 
that shape effect deviations are not effectively removed by the 
three averaging devices we have introduced in Section II C, 
but fortunately these deviations are very small. 



C. Twist Boundary Conditions Averaging 

For a system of interacting spinless fermions, we cannot di- 
rectly compute the exact infinite-system cluster DM. We can, 
however, choose to work with (i) an approximate ground-state 
wave function of an infinite system, or (ii) the exact ground- 
state wave function of a finite system, or (iii) an approximate 
ground-state wave function of a finite system. As reported in 
Section IV, we chose option (ii), where the exact ground-state 
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FIG. 3: Zero- and one-particle weights of the cluster DM of a 5-site, 
cross-shaped cluster for systems of noninteracting spinless fermions 
with periodic boundary conditions imposed. The finite (4, - l)x(l, 4) 
(•) and (8, -2) x (2, 8) (■) systems have the same shape but different 
sizes. The solid curves are the zero- and one-particle weights for the 
infinite system. For w si , we see the finite domain effect deviations for 
the (4, -l)x(l,4) system is practically gone by the time we get to the 
(8, -2) x (2, 8) system. For the rest of the one-particle weights, the 
finite domain effect is largely removed in the (8, -2) x (2, 8) system, 
but shell effect persists. In fact, shell effect in Wd is still visible when 
we go to the (32, -8) x (8, 32) system (not shown), which has N = 
1088 sites. 



wave function is obtained through numerical exact diagonal- 
ization. 

Exact diagonalization severely limits the sizes of the finite 
systems we can work with (see Appendix of Ref.9 for formula 
on size of Hilbert space for the strongly-interacting model 
given by (1.1)). With aggressive memory reduction measures, 
it is possible (but not necessarily feasible) to exactly diagonl- 
ize finite systems with up to 30 sites for all filling fractions. 
However, as we have illustrated in Section III B, for systems 
so small, the numerical cluster DM spectra would be plagued 
by strong finite size effects, most notably by the shell effect. 
This is where twist boundary conditions averaging comes in. 

Very crudely speaking, we can think of averaging numeri- 
cal observables over M twist boundary conditions for a finite 
system with N sites as being equivalent to computing these 
numerical observables for a single finite system of N* > N 
sites, subject to only periodic boundary conditions. In the 
best-case scenario, the effective system size N* might be as 



8 



large as MN, though N* will typically grow slower than O(M). 
The computational cost of performing exact diagonalization 
for a system with N sites over M twist boundary conditions 
is on the order of 0(MN 3 ), whereas the computational cost 
of performing exact diagonalization just once for a system 
with N* > N sites is 0(N* 3 ). So long as the effective system 
size N* grows faster than 0(M 1 ^), it would be computational 
cheaper to employ the method of twist boundary conditions 
averaging, instead of exactly diagonalizing a single large finite 
system, to reduce finite size effects. The detail dependence of 
N* on M will of course depend on the nature of the observable 
of interest. 

From the detailed study undertaken in Appendix D of 
Ref. 1 1, we know that there are cuts and cusps on the twist sur- 
face {¥(<px, 0->,)l < 2| v P(0*> <Pv)) °f a generic observable O, where 
\¥(<p x , <py)) is the many-body ground state of a finite A^-site 
system subject to twist boundary conditions with twist vector 
= ((p x , 4>y). For non-square systems, these cusps and cuts 
demarcate features with a hierarchy of sizes on the twist sur- 
face. The 'typical' twist surface feature has a linear dimension 
of 2n/ y/N. These are decorated by fine structures with lin- 
ear dimension 2n/N, which are in turn decorated by hyperfine 
structures with linear dimensions In/N 2 . The number of in- 
tegration points we must use is therefore determined by what 
feature size we want to integrate faithfully. 

For the purpose of this study, we decided to integrate the 
fine structure on the twist surface faithfully. Therefore, we 
chose to average the spectrum of pc over a q — 16 Monkhorst- 
Pack grid (which consists of 256 integration points in the First 
Brillouin Zone) for the (4, -1) x (1,4) system with N - 11 
sites. We find that twist boundary conditions averaging does 
indeed result in an averaged spectrum which approximates the 
infinite-system limit well (see Fig. 4). This averaging device, 
however, does not completely eliminate shell effects, as can 
be seen from the twist boundary conditions-averaged Wd(n). 
To reduce the bias this creates for one particular choice of 
finite system, we combined the twist boundary conditions- 
averaged spectrum of pc for various finite systems. This 
is shown in Fig. 5. We will overlay the cluster DM spectra 
from several finite systems in the same way, to derive a twist 
boundary conditions-averaged approximation to the infinite- 
system cluster DM spectrum of strongly-interacting spinless 
fermions. 



IV. STRONGLY-INTERACTING SYSTEM 

In this section, we compute the cluster DM for interacting 
spinless fermions. As with the case of noninteracting spinless 
fermions, the cluster DM evaluated directly from the ED of 
various finite systems are severely affected by finite size ef- 
fects. In Section III, we saw how finite size effects can be 
significantly reduced when the cluster DM is averaged over 
various twist boundary conditions, in the sense that the av- 
eraged cluster DM weights from different finite systems at 
various filling fractions fall close to their respective infinite- 
system limits. Applying twist boundary conditions averaging 




FIG. 4: One-particle weights of the cluster DM of a 5-site, cross- 
shaped cluster within systems of noninteracting spinless fermions. 
The performance of twist boundary conditions averaging, using q = 
16 Monkhorst-Pack special-point integration, in reducing finite size 
effects for the (4,-1) x (1,4) system (•) is checked against the 
(4,-1) x (1,4) (O) system with periodic boundary conditions im- 
posed. The solid curves are the cluster DM weights from the infinite 
system. 



onto the interacting cluster DM, we find that finite size effects 
are reduced, but not as dramatically as for the noninteracting 
cluster DM. Nevertheless, the averaged cluster DM weights 
from different finite systems with different number of particles 
are sufficiently consistent with each other that we can plot a 
smooth curve interpolating the averaged cluster DM weights. 

As explained in Section I, our interest in studying the 
strongly-interacting model (1.1) of spinless fermions with in- 
finite nearest-neighbor repulsion is to understand how the 
cluster DM evolves with filling, given that we expect in 
this model crossovers between regimes of qualitatively differ- 
ent states. Furthermore, we had proposed an operator-based 
method of truncation which was justified by the fact that the 
noninteracting cluster DM is generated from a set of single- 
particle operators. Since we proposed the use of this trunca- 
tion scheme for interacting systems, we are interested to know 
whether the structure of the interacting cluster DM is such that 
it can also be generated, perhaps approximately, from a set of 
single-particle operators. 

To this end, we present in Sections IV A, IV B and IV C, 
the zero-, one- and two-particle cluster DM weights of the 
strongly-interacting system of spinless fermions. We check 
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about the correlations in the strongly-interacting ground state, 
and so we move on to consider the one-particle cluster DM 
weights. 



FIG. 5: One-particle weights of the cluster DM of a 5-site, cross- 
shaped cluster, for the (3, -2)x (2, 3) (•), (4, l)x (1, 3) (O), (4, -l)x 
(1, 3) (□), (4, -l)x(l,4) (■), (4, -l)x(2, 3) (0), (4, -2)x(2,4) (♦) and 
(5, 1) x (1, 5) (A) systems of noninteracting spinless fermions subject 
to twist boundary conditions averaging, using q = 16 Monkhorst- 
Pack special-point integration. Also shown are the cluster DM 
weights of the infinite system. 




FIG. 6: Zero-particle weight of the cluster DM of a 5-site, cross- 
shaped cluster, for the (3,-2) x (2,3) (•), (4,1) x (1,3) (O), 
(4, -l)x(l,3) (□), (4, - l)x(l,4) (■), (4,-l)x(2, 3) (0) and (4,-2)x 
(2, 4) (♦) systems of strongly-interacting spinless fermions subject to 
twist boundary conditions averaging, using q = 8 Monkhorst-Pack 
special-point integration. At n = and n = \, we know analytically 
that wo = 1 and w = respectively, and the solid 'curve' interpo- 
lates between these two known limits and the equally weighted data 
points at finite filling fractions < n < \. Also shown as the dashed 
curve is the zero-particle weight of the infinite system of noninter- 
acting spinless fermions. 



whether it is possible to: (i) write the two-particle eigen- 
states as the product of one-particle eigenstates; and (ii) pre- 
dict the relative ordering of the two-particle weights based on 
the relative ordering of the one-particle weights. We then dis- 
cuss in Section IV D whether these two criteria are met, and 
whether it is feasible to design an operator-based DM trun- 
cation scheme, similar to the one described in Ref.7, for the 
strongly-interacting system. 

A. Zero-Particle Cluster DM Weight 

The zero-particle cluster DM weight calculated for various 
finite strongly-interacting systems is shown in Fig. 6. Also 
shown in the figure is the zero-particle cluster DM weight of 
the infinite system of noninteracting spinless fermions. As we 
can see, the zero-particle weights of the respective systems 
only start differing significantly from each other for « > 0.1. 
With repulsive interacting between spinless fermions, it is 
more difficult in a congested system (n > 0.2) to form an 
empty cluster of sites from quantum fluctuations. As a re- 
sult, the strongly-interacting Wq falls below the noninteracting 
wq. However, this fact alone does not tell us anything more 



B. One-Particle Cluster DM Eigenstates and Their Weights 

In our Operator-Based DM Truncation Scheme developed 
in Ref.7 for a noninteracting system, the one-particle clus- 
ter DM weights play a very important role, since we select 
which one-particle operators to keep or discard based on the 
negative logarithm of these numbers. We expect the one- 
particle cluster DM weights, though not entirely sufficient by 
themselves, would also play an important role in defining an 
operator-based truncation scheme. Therefore, in this section, 
we present results for a series of calculations to determine the 
infinite-system limit of the one-particle cluster DM spectra for 
our strongly-interacting system as a function of filling fraction 
h, and discuss their implications for an operator-based DM 
truncation scheme. 

Though we really do need to worry about the evolution of 
the structure of and \s<i) as a function of « in both the non- 
interacting and strongly-interacting systems, the one-particle 
weights are ordered by their magnitudes as w Sl > w p > 
Wd > w Sl for both systems. But while the noninteracting one- 
particle weights go down by roughly one order of magnitude 
as we go through the sequence w S] — > w p — > Wj — > w Sl , we 
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FIG. 7: One-particle weights of the cluster DM of a 5-site, cross- 
shaped cluster within a system of strongly-interacting spinless 
fermions. 



see from Fig. 7 that the interacting one-particle weights decay 
more slowly along this same sequence. 

We studied the finite (3,-2) x (2,3) (•), (4,1) x (1,3) 
(O), (4,-1) x (1,3) (□), (4,-1) x (1,4) (■), (4,-1) x (2,3) 
(0) and (4, -2) x (2,4) (♦) systems subject to twist boundary 
conditions averaging, using q — 8 Monkhorst-Pack special- 
point integration. At a filling fraction of n = 0, the sys- 
tem approaches the noninteracting limit, and thus all the one- 
particle weights are zero. At half-filling n = \, the two- 
fold degenerate checker-board ground state is unaffected by 
twist boundary conditions averaging. We can thus perform 
degeneracy averaging analytically, to find that w Sl = I and 
w p - Wd - w S2 = 0. The solid 'curves' in Figure 7 interpo- 
late between these two known limits and the equally weighted 
data points at finite filling fractions < n < A. Also shown in 
Figure 7 as the dashed curves are the one-particle weights of 
the infinite system of noninteracting spinless fermions. 

When only periodic boundary conditions are imposed, there 
is significantly more 'scatter' in the one-particle weights as a 
function of filling fraction h, for interacting systems of dif- 
ferent sizes, than for noninteracting systems of different sizes. 
Averaging the one-particle weights of the interacting systems 
over various twist boundary conditions visibly reduces this 
'scatter', even though the remnant 'scatter' seen in Figure 7 is 
still rather large, compared to case for the noninteracting sys- 
tem (Figure 5). From our own detailed study of the method 



of twist boundary conditions averaging for noninteracting 
systems, 11 we know that twist boundary conditions averaging 
effectively removes finite size effects from some observables, 
but not for others. We have no reason to expect twist bound- 
ary conditions averaging to apply equally effectively over the 
same set of observables, when we go from noninteracting 
systems to interacting systems. Conversely, observables for 
which twist boundary conditions averaging is ineffective in 
noninteracting systems, might be effectively twist-boundary- 
conditions-averaged in interacting systems. With only input 
from the exact diagonalization of finite systems, and with- 
out employing system-size extrapolations, the method of twist 
boundary conditions averaging offers us the best hope of gain- 
ing insight into the infinite-system properties we seek. 

We expect that the remnant 'scatter' in the twist-boundary- 
conditions-averaged one-particle weights will be reduced, if 
we had not made the nearest-neighbor repulsion infinite. 
There are two reasons why we did not also study the case 
of t <k V < co. First, for a fixed system size N and parti- 
cle number P, the Hilbert space for the V < oo system would 
be much larger than that of the V — > oo system. A parallel 
study for V < oo of the V — > oo system sizes and particle 
numbers reported in this paper, with twist boundary condi- 
tions averaging, will require unacceptably long total compu- 
tation time. Second, and more importantly, we believe that as 
long as V/t < oo is large, the qualitative implications on the 
Operator-Based DM Truncation Scheme would be similar to 
the case where V/t — > oo, and therefore it suffices to examine 
the latter case, which is computationally much more manage- 
able. In any case, we do not believe the remnant 'scatter' 
in the twist-boundary-conditions-averaged DM weights will 
hamper our efforts in drawing qualitative conclusions regard- 
ing the applicability, or otherwise, of the Operator-Based DM 
Truncation Scheme for interacting Fermi liquids. 

In the Operator-Based DM Truncation Scheme described 
in Ref.7, we discard one-particle cluster DM eigenstates with 
very small weights, and keep only the many-particle cluster 
DM eigenstates built from the retained one-particle eigen- 
states. The sum of weights of the truncated set of cluster DM 
eigenstates will then be very nearly one, if the discarded one- 
particle weights are all very small compared to the maximum 
one-particle weight. As we can see from Fig. 7, the ratio of the 
largest one-particle weight, w Sl , to the smallest one-particle 
weight, w Sl , is not large enough for us to justify keeping 
and discarding \s2), except when the system is very close to 
half-filled. 

We believe that the one-particle cluster DM weights are so 
close to each other in magnitude, because of the net 'strength' 
of interactions straddling the cluster and its environment is 
strong compared to the net 'strength' of interactions strictly 
within the cluster. Unfortunately for our five-site cluster, 
which was chosen because it is the smallest non-trivial cluster 
having the full point group symmetry of the underlying square 
lattice, 16 the sites within the cluster are poorly connected, i.e. 
a cluster site is on average connected to more environmental 
sites than to other cluster sites. To have more of the total inter- 
actions of the cluster with the system be confined within the 
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cluster, a cluster significantly larger than the five-site cluster 
studied in this paper will be needed. This large cluster must 
then be embedded in a finite system that is larger still, making 
exact diagonalization studies unfeasible. 



C. Two-Particle Cluster DM Eigenstates and Their Weights 




FIG. 8: Two-particle weights of the cluster DM of a 5-site, 
cross-shaped cluster within a system of strongly-interacting spinless 
fermions. 

While it is desirable to have a broader distribution of one- 
particle weights, our more important task is to examine how 
closely the many-particle cluster DM eigenstates can be ap- 
proximated as products of one-particle cluster DM eigen- 
states. In particular, we look at the two-particle cluster DM 
eigenstates, and find that of the two-particle states listed in 
(2.27), the only states which are allowed by the no-nearest- 
neighbor constraint to appear in the cluster Hilbert space are 
\S), \P±), \P'i) and \D). We know therefore that the two- 
particle sector of pc comprises a 1 X 1 5 -diagonal block (with 
weight ws ), a 1 x 1 D-diagonal block (with weight wd), and 
two degenerate 2x2 P-diagonal blocks (with weights vvp, and 
wp n ). The two-particle weights are shown as a function of 
filling h in Fig. 8. 

For the finite (3,-2) x (2,3) (•), (4,1) x (1,3) (O), 
(4, -1) x (1, 3) (□), (4, -1) x (1,4) (■), (4, -1) x (2, 3) (0) and 
(4, -2) x (2, 4) (♦) systems studied, subject to twist boundary 
conditions averaging, using q — 8 Monkhorst-Pack special- 



point integration, all the two-particle weights are zero at h = 
as the systems approach the noninteracting limit. At half- 
filling « = 5, we again perform degeneracy averaging analyt- 
ically on the two-fold degenerate checker-board ground state 
to find that all the two-particle weights are zero. In Figure 8, 
the solid 'curves' interpolates between these two known limits 
and the equally weighted data points at finite filling fractions 
< n < i. 

We realized that there are significantly fewer nontrivial two- 
particle eigenstates of pc than predicted from the combination 
of one-particle eigenstates. From the point of view of im- 
plementing the Operator-Based DM Truncation Scheme, this 
poses no problem ;/ the non-occurring two-particle states are 
predicted to have small enough weights that they will be ex- 
cluded by the truncation scheme. However, we find that this 
is not the case. For example, the two-particle state \S'), which 
does not occur, is predicted by simple combination of the one- 
particle states \s) and \s) to have a weight comparable to that 
of the two-particle state \S), which does occur. 

Of the two-particle weights that are allowed by the no- 
nearest-neighbor constraint, we expect their weights to fol- 
low the sequence wp t > wd > ws, if they can indeed to 
thought of products of one-particle states. From Fig. 8, we 
indeed observe this sequence of two-particle weights, even 
though their actual magnitudes (calculated as the product of 
one-particle weights divided by the zero-particle weight) do 
not come out right. This observation is encouraging, because 
we might yet be able to push a variant of the Operator-Based 
DM Truncation Scheme through, by introducing constraints 
on how many-particle cluster DM states can be built up from 
one-particle cluster DM states. 

D. Signatures of Fermi-Liquid Behaviour in the Cluster DM 

Over broad ranges of filling fractions, the ground state of 
our strongly-interacting model (1.1) of spinless fermions is 
expected to be an interacting Fermi liquid. While we under- 
stand the cluster DM structure of a noninteracting Fermi liq- 
uid completely, we do not yet understand how an interacting 
Fermi liquid will manifest itself in the structure of its clus- 
ter DM. Unlike a noninteracting Fermi liquid, the ground- 
state DM of an interacting Fermi liquid will not simply be 
the exponential of a noninteracting pseudo-Hamiltonian. Nev- 
ertheless, we expect that the interacting pseudo-Hamiltonian 
H appearing in the interacting Fermi-liquid ground-state DM 
p = exp(H) can be made to look like the sum of a noninteract- 
ing pseudo-Hamiltonian Hq, and a much weaker interaction 
term Hi , by a canonical transformation. From Landau's Fermi 
liquid theory, we know that such a canonical transformation 
(similar to the one which relates Landau quasi-particles to 
bare fermions) works by burying much of the bare interac- 
tions within the quasi-particles. 

In tracing down the ground-state DM, our hope then is that 
the cluster DM can also be written, after a canonical transfor- 
mation local to the cluster, as the exponential of the sum of 
a noninteracting pseudo-Hamiltonian Hc,o (which is perhaps 
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related to Hq in the same way as for the noninteracting Fermi 
liquid), and a weak interaction term Hc.i- We suspect that the 
criterion for this to be possible is that we must be able to con- 
struct approximate quasi-particles using only cluster states to 
absorb the bare interactions. However, we also believe that 
the requirement that the canonical transformation act strictly 
within the cluster will fail to completely incorporate interac- 
tions that straddle the cluster and its environment. 

In Section IV C, we found that for our strongly-interacting 
system, the two-particle cluster DM eigenstates look nothing 
like simple products of two one-particle cluster DM eigen- 
states each. In fact, many combinations of two one-particle 
cluster DM eigenstates give invalid two-particle cluster states 
that violate the condition of no-nearest-neighbor occupation. 
It is tempting, based on this observation, to say then that the 
cluster DM is not the exponential of an approximately nonin- 
teracting pseudo-Hamiltonian. However, it must be remem- 
bered that in an interacting Fermi liquid, the quasi-particles 
are also not single bare particles. Instead, they are superpo- 
sitions of states containing different number of bare particles, 
which leads us to think of a quasi-particle as a bare particle 
being screened by other bare particles in its vicinity. 

With this in mind, we realized that to identify the quasi- 
particle structure of the cluster DM, we need to construct 
appropriate linear combinations of the P-particle cluster DM 
eigenstates, so that the cluster DM, when written in terms of 
these 'quasi-particles', look like the exponential of a noninter- 
acting pseudo-Hamiltonian Hc.o plus a weak interaction term 
Hc,\- This involves writing the pseudo-Hamiltonian He as 
a sum of terms, representing the independent quantum fluc- 
tuations associated with each of the quasiparticles. This can 
be accomplished by defining an operator singular value de- 
composition of the cluster DM with respect to an appropriate 
operator norm, which forms the basis forjudging whether the 
quantum fluctuations associated with two linear combination 
of bare operators are independent. Details of such an opera- 
tor singular value decomposition will be reported in a future 
paper. 15 

V. SUMMARY & DISCUSSIONS 

To summarize, we have calculated numerically the clus- 
ter DM for a cross-shaped cluster of five sites within both a 
system of noninteracting spinless fermions described by the 
Hamiltonian (3.1), and a system of strongly-interacting spin- 
less fermions described by the Hamiltonian (1.1). For the 
noninteracting system, the cluster DM was obtained from the 
cluster Green-function matrix using the exact formula (3.2) 
obtained in Ref.5, whereas for the interacting system, the clus- 
ter DM was obtained from the exact-diagonalization ground- 
state wave function by tracing down degrees of freedom out- 
side of the cluster. For the purpose of making the comparison 
of the cluster DM spectra more straightforward, we worked 
with the same collection of finite non-square systems for in- 
teracting and noninteracting spinless fermions. 

To make the results of this comparison less dependent on 



the geometry of the finite systems chosen, degeneracy averag- 
ing followed by orientation averaging of the cluster DM spec- 
tra were carried out. When combined, these two averaging 
apparatus has the effect of restoring full square-lattice sym- 
metry to the cluster DM. We also analyzed in detail the finite 
size effects not removed by degeneracy and orientation aver- 
aging, by inspecting the numerical cluster DM spectra coming 
from finite noninteracting systems of various sizes. By in- 
creasing the system size systematically, we find visually that 
the infinite-system limit of the cluster DM is 'attained' when 
the system size reaches a couple of hundred sites. 

Noting that these are forbiddingly large system sizes to 
work with for interacting systems, where the ground-state 
wave functions have to be obtained via exact diagonaliza- 
tion, we then tested the apparatus of twist boundary condi- 
tions averaging on finite systems with between 10 and 20 
sites. For noninteracting spinless fermions, we find that the 
twist boundary conditions-averaged cluster DM weights for 
different finite systems and different filling fractions indeed 
fall close to the various infinite-system limits. Since we do not 
perform system-size extrapolations, we interpolate between 
the degeneracy-, orientation-, and twist boundary conditions- 
averaged cluster DM weights for the various finite interacting 
systems and their respective accessible filling fractions, and 
take the result curves to be our best approximation of the clus- 
ter DM spectrum of the infinite interacting system. 

Comparing the twist boundary conditions-averaged cluster 
DM spectra for the noninteracting and strongly-interacting 
systems, we find similar qualitative behavior in the zero- 
particle weights as functions of filling fraction, and quali- 
tatively different behaviours in the one-particle weights as 
functions of filling fraction. However, the relative ordering 
w Sl > w p > Wd > w s , is the same at all h < | for both 
systems. Quantitatively, we find for noninteracting spinless 
fermions that the one-particle weights go down by roughly 
one order of magnitude each time as we go through the se- 
quence w Sl — > w p — > Wd — > h>, 2 . For strongly-interacting 
spinless fermions, the one-particle weights decay much more 
slowly along the sequence. 

The implications this observation have for the Operator- 
Base DM Truncation Scheme developed in Ref.7 is that, for a 
small fixed fraction of one-particle eigenstates retained, the to- 
tal cluster DM weight of eigenstates retained would be much 
smaller for the strongly-interacting system compared to the 
noninteracting system, since the ratio of the smallest to the 
largest one-particle weights, w Sl /w Sl , is not very much smaller 
than one. This narrow distribution of one-particle cluster 
DM weights aside, we observed that the relative ordering 
Wp, > wd > ws of the two-particle cluster DM weights, pre- 
dicted based on the combination of one-particle cluster DM 
weights, is confirmed numerically — even though the pre- 
dicted weights are off. This suggests that we might be able 
to push a variant of the naive Operator-Based DM Truncation 
Scheme through, by introducing additional rules on how one- 
particle cluster DM eigenstates can be combined to give only 
the valid many-particle cluster DM eigenstates. We did not 
attempt to implement such a truncation scheme, and test how 
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badly its numerical accuracy is affected by the ratio w Sl /w Sl 
(by calculating the dispersion relation, for example, as was 
done in Ref.7), because we feel that such a naive scheme was 
not in the spirit of finding an appropriate 'quasi-particle' de- 
scription for the low-energy excitations of our system of inter- 
acting spinless fermions given by (1.1). 

Finally, we realized that our numerical studies do not al- 
low us to conclude whether the cluster DM of the interact- 
ing system furnishes a good 'quasi-particle' description for 
the strongly-interacting system. To be able to check this, we 
must be able to construct appropriate superpositions of clus- 
ter DM eigenstates with different particle numbers, so that 
the pseudo-Hamiltonian He ~ log pc looks like the sum of a 
noninteracting pseudo-Hamiltonian i?c,o and a weak interact- 
ing pseudo-Hamiltonian Hc,i- Instead of a simple eigenvalue 
problem for the cluster DM, the problem of discovering what 
'quasi-particle' operators make up the cluster DM is an oper- 
ator singular value decomposition problem. We will carefully 
define this operator singular value decomposition in a future 
paper, and describe how it can be applied to the density ma- 
trix of two disjoint clusters to systematically extract operators 
associated with independent quantum fluctuations within each 
cluster, and their inter-cluster correlations. 15 
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APPENDIX A: COMPUTATIONAL COMPLEXITY OF 
CLUSTER DENSITY-MATRIX CALCULATION 

In this appendix we compare the naive algorithm and the 
presorted inner product algorithm, based on (2.9) and (2.10) 
respectively, for numerically computing the cluster density 
matrix, and determine their computational complexities. To 
begin, we denote by D(P) the size of the system Hilbert space 
with P particles, Dc(Pc) the size of the cluster Hilbert space 
with Pc particles, and De(Pc) the size of the environment 
Hilbert space with Pe = P - Pc particles. Noting that there 
can be no matrix elements between cluster configurations with 
different number of particles, we calculate each Pc sector of 
the cluster DM separately. To keep our notations compact, 
let us drop the P and Pc dependences in D(P), and Dqe(Pc) 
respectively from this point onwards, and reinstate these de- 
pendences only when necessary. Readers are referred to Ap- 
pendix A in Ref. 1 1 for more technical details on the compu- 
tational implementation of this trace-down calculation of the 
cluster DM. 



In the naive algorithm based on (2.9), the cluster DM ma- 
trix elements are computed by starting nested 'for' loops in 1 
and 1', each running over Dc indices. For each pair of clus- 
ter configurations |1) and |1'), one would need to then deter- 
mine which of the P-particle configurations |j) contain the 
two cluster configurations. This involves running through 
the D configurations in the system Hilbert space, and for 
each configuration, comparing the P occupied sites with the 
Pc occupied sites in the cluster configurations |1) and 
The computational effort incurred for this matching is thus 
O(DP). Two vectors of indices, i, whose entries are the in- 
dices of system configurations |j) giving cluster configuration 
and i', whose entries are the indices of system configu- 
rations j) giving cluster configuration |1'), are obtained. The 
lengths of these index vectors vary, but are of O(De). One 
can then compare the two index vectors, at a computational 
cost of 0(D 2 E ), to find which pairs of system configurations 
giving cluster configurations |1) and |1') share the same envi- 
ronment configuration. Following this, one can sum over the 
amplitude of such pairs, at a computational cost of O(De), 
to obtain the cluster DM matrix element (l|pc|l'). For this 
naive algorithm, the net computational effort is on the order 
of D 2 (DP + D 2 E + D E ) ~ D 2 C (DP + D 2 E ). 

In the pre-sorted inner-product algorithm based on (2.10), 
we need to first run through D system configurations to pre- 
sort the amplitudes in the ground-state wave function. For 
each system configuration |j), we determine at a computa- 
tional cost of P comparisons what cluster configuration |1) and 
environment configuration |m) it contains. We then search 
through the cluster and environment Hilbert spaces to deter- 
mine what the indices of |1) and |m) are in their respective 
Hilbert spaces, which incurs a computational effort on the 
order of DcPc and DePe respectively. Once these indices 
are determined, the amplitudes in the ground-state wave func- 
tion are organized into a Dc x De matrix. The net compu- 
tational expenditure is thus on the order of D(P + DcPc + 
D e Pe) ~ D(D c Pc + D e Pe)- After sorting the ground-state 
wave function, we can then start nested for loops in 1 and 1', 
each running over Dc indices, to evaluate the matrix element 
(l|pc|l') as the inner product between two vectors of length 
De- This trace-down stage incurs a computational cost of 
0{D 2 c De)- Overall, the computational cost is on the order of 
D(D C P C + D E P E ) + D 2 C D E . 

For models allowing nearest-neighbor occupation, the sys- 
tem Hilbert space is the direct product of the cluster Hilbert 
space and the environment Hilbert space, i.e. D = DcDe- 
Since the number P of particles is small in any reasonable 
exact diagonalization, we can treat it as a 0(1) constant. For 
small clusters, the size Dc of the cluster Hilbert space will 
also be small, so that the size De of the environment Hilbert 
space will be comparable in magnitude to D. With these con- 
siderations, we find that the computational cost for the naive 
algorithm is 0(DD 2 c + D 2 C D 2 E ) ~ 0(D 2 ), while the computa- 
tional cost for the inner-product algorithm with pre-sorting is 
0(D + DD E + D E ) ~ 0(D 2 ). The efficiency of the two al- 
gorithms therefore depend on the prefactors, the estimation of 
which requires more thorough analyses of the two algorithms. 
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For a model such as (1.1), where nearest-neighbor occupa- 
tion is forbidden, the system Hilbert space is smaller than the 
direct product of the cluster Hilbert space and the environment 
Hilbert space, i.e. D < DqDe- Given again that P and Dq are 
small numbers, the computational cost for the two algorithms 
are essentially determined by the ratio DcDe/D. This ratio 
is strongly dependent on the dimensionality of the problem: 
the superfluous configurations generated by the direct prod- 
uct of the cluster Hilbert space and the environment Hilbert 
space are invalid because they contain nearest-neighbor sites, 
right at the boundary between the cluster and its environment, 
which are occupied. In one dimension, the number of super- 
fluous configurations is small, because the boundary between 
the cluster and its environment consists only of two bonds, 



whatever the size Nc of the cluster. In two dimensions, the 
boundary between the cluster and its environment is a line 
cutting roughly V^Vc bonds. The number of superfluous con- 
figurations is then proportional to exp(a2 V^c)> where 0-2 is 
a constant prefactor which depends on the shape of the clus- 
ter. In d dimensions, the number of boundary bonds is on the 
order of N^,~ , and the number of superfluous states is pro- 
portional to exp(arrf./V^ ). Therefore, in dimensions greater 
than one, DqDe become increasingly larger than D as Nc is 
increased, and the inner-product algorithm with pre-sorting, 
which involves only one power of DqDe, is more efficient 
than the naive algorithm, which involves (DcAe) 2 . 
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